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We discuss the origin and relevance for computer simulations of a strong finite-size effect that 
appears when using the Ewald summation formula. It can be understood as arising from a volume- 
dependent shift of the potential in a finite, periodic box relative to the infinite volume limit. This 
shift is due to the fact that the "zero of energy" for a periodic system cannot be defined by letting 
the interacting particles be separated by an infinite distance; the correct definition corresponds to 
■ setting its k = Fourier mode to zero. The implications of this effect for computer simulations are 

discussed. 



in 



(N 
> 



in 



I. INTRODUCTION 



The correct and efficient treatment of the long-range 
electrostatic interactions in simulations of aqueous so- 
lutions continues to be an issue of intense research. 
The Ewald summation technique ||-|3| appears to be a 
promising approach for the consistent treatment of these 
• interactions. Several groups have discussed that the trun- 
cation of the electrostatic forces, either spherical or using 
the "minimum image" convention, results in severe mod- 
ifications of the dielectric behavior There are in- 
dications that these modifications have quite noticeable 
effects on the behavior of aqueous solutions in computer 
simulations P-p"3|. 

In spite of this, the Ewald summation method has 
not been extensively applied to simulations of macro- 
molecules in solution, where the preferred method con- 
tinues to be spherical truncation of the Coulomb poten- 
tial. With the rapid advances in computer technology 
and new, improved algorithms that speed the calcula- 
. , tion of the long-range forces without sacrificing accuracy 
0-|l8|, this situation may soon change. Thus, a discus- 
sion of the artifacts that affect the computer simulations 
is timely. In this paper we focus our attention on a finite- 
size effect, that is, an artifact that arises from the fact 
that the system being simulated is constrained to a finite 
volume, V. This effect can be quite strong and, although 
it disappears in the limit V — ► oo, the convergence to this 
limit is not particularly rapid. It should be noted, how- 
ever, that volume-dependencies affect any simulation of 
a finite system, and are not just an artifact of using the 
Ewald summation method, although the Ewald method 
lends itself to a more systematic theoretical treatment. 

The effect described in this paper can be seen from 
two different perspectives: either as the result of the fact 
that constraining the system to a finite box introduces 
a volume-dependent natural scale for the interaction en- 
ergy of a system of two ions; or as the dependence on 
the unit cell volume, V, of the electrostatic interaction 
energy of a pair of charges immersed in a continuum of 
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dielectric e (with periodic boundary c ond itions). These 



two views are expounded in sections III and IV. In a 
related paper Del Buono et al. demonstrate by ex- 
plicit computer simulations in water that this effect is 
not negligible. 

The structure of the paper is as follows. In the second 
section we discuss the most general form of the electric 
potential in a periodic box of volume V and show that 
simple physical requirements lead to the (tinfoil) Ewald 
potential augmented by an "extrinsic" potential, propor- 
tional to the unit cell's dipole moment. This section does 
not contain new results but can shed new light on old re- 
sults that are still not quite clear in the technical litera- 
ture |2(J . In section III the classical electrostatics expres- 
sion for the interaction energy of two point charges in a 
periodic box is discussed. It is shown that the interaction 
energy coincides with the potential only if the latter satis- 
fies the condition that its volume average vanishes. This 
condition generalizes the usual "zero potential at infin- 
ity" rule in infinite systems. Section IV presents results 
of molecular dynamics simulations of two ions in aqueous 
solution in order to illustrate the issue of volume depen- 
dence and its importance to computer simulations. We 
focus on the free energy of charging a system of two ions 
in a bath and show that, again, the condition of vanishing 
volume average arises quite naturally. Next we discuss a 
numerical example to put into perspective the size of the 
effect. Finally the relevance to computer simulations of 
aqueous solutions is discussed. 



II. THE LATTICE SUMS 

Since the seminal work of de Leeuw et al. |^l|,^2| it has 
been known that the lattice summation method, that is, 
the summation of the potentials and forces produced by 
an infinite number of replicas of a system of charges (unit 
cell), produces not one but a whole family of potentials, 
which differ in a term proportional to the dipole moment 
d of the unit cell and take the form 
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Although for A ^ this potential is not periodic due 
to the term linear in a;, it nevertheless gives rise to a pe- 
riodic electric field, which is the only constraint that can 
be imposed on physical grounds. The constant of pro- 
portionality, A, can be related to the dielectric constant 
£rf of a medium surrounding a large sphere of repli- 
cas of the unit cell. In the limit of an infinite sphere a 
residual dependence on this dielectric constant remains, 
giving rise to this "extrinsic" potential |23|| . Sometimes 
the choice of £rf is called the "boundary condition," not 
to be confused with the fact that these models impose 
periodic boundary conditions on the particles. The value 
£rf = oo (tinfoil or conducting boundary) is very special 
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since it alone results in a potential that is periodic due 
to the absence of this dipolar term. All other choices 
necessarily break periodicity and their use in molecular 
dynamics simulations introduces some subtle questions. 
In what follows when we talk about the "Ewald poten- 
tial" we mean that which results from the choice of tin- 
foil boundary conditions. If other "boundary conditions" 
are chosen another finite-size effect related to this dipolar 
term appears. This will not be discussed further in this 
paper. 

In their excellent mathematical treatment de Leeuw 
et al. |^l],^H were concerned with the calculation of the 
conditionally convergent (i.e., ill-defined) sum 



\xi — Xj — Ln\ 



that defines the potential energy of a system of charges 
in a periodic box of side length L (notice that for n = 
the diagonal terms i — j are dropped). Formally this 
sum can be rewritten as 

l^2liQjH x i ~ x j) ( 3 ) 
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where 



n " 11 

It is this sum, or rather a regularization of it, that we 
will be concerned with. 

First of all, note that if the system is electrically neu- 
tral we can add to $ any constant term without changing 
the potential energy since 



1,3 




QiQj x c = c x qA = 0. (5) 



The second observation is that §(xi — Xj) diverges as 
1/fiEj — Xj\\ when Xi — > Xj. This short-distance diver- 
gence is well known and has nothing to do with the fact 
that the sum in equation (||) is ill-defined; it is just an 
artifact of using point charges, which have a divergent 
self-energy even in an infinite space. Thus the diagonal 
sum i = j should be redefined to subtract this trivial di- 
vergence. In any case, the diagonal sum contributes only 
to the self-energy of the system and does not influence 
the dynamics. 

The third observation is the following: the sum in 
equation (||) diverges. In their treatment de Leeuw et 
al. actually found that the divergence is independent of 
x, and there is no effect for an electrically neutral sys- 
tem. Once this divergence is removed what remains is 
conditionally convergent. We should mention that there 
seems to be a confusion about this point: the only re- 
quirement that needs to be made to obtain sensible re- 
sults is that the system be electrically neutral and not, 
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as has been mentioned in the literature [£(J , that the net 
dipole moment vanish as well. If there is a net dipole mo- 
ment the sum in equation (||) is conditionally convergent 
(not divergent), which means that there is not a unique 
answer. This nonuniqueness reflects the fact that some 
properties of a system depend not only on the intrinsic 
characteristics of the system but also on the shape and 
other boundary conditions. Due to its long-range na- 
ture, the electrostatic energy is a prime example of this 
dependence. A clear discussion of this point, albeit quite 
technical, can be found in f p3[ . 

Since $(a; — u) is the electric potential produced by a 
unit charge at u (plus all its infinite copies), we might 
start by requiring that $ satisfies Poisson's equation |2J] 

V 2 $(cc -u) = -4tt5 p (x - u) (6) 

and that $ be a periodic function: $(a; + Ln) = $>(x) 
(here S p is the periodic Dirac 8- function). This is, how- 
ever, not consistent, as integrating both sides over the 
unit cell shows. The left hand side gives identically zero, 
as can be seen by using Stoke's theorem to convert the 
volume integral to a surface integral and the fact that 
the electric field is periodic, while the right hand side 
does not vanish. This is simply a consequence of the 
well known fact that in a periodic system the net charge 
must vanish. We can easily satisfy the neutrality require- 
ment by adding a homogeneous background charge den- 
sity that exactly cancels the unit charge in equation (0): 



V 2 $(a; - u) = -Att ( 5 p (x - u 
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where we have also written the right hand side in Fourier 
space (here V = L 3 ). 

With hindsight we can find a particular solution of 
equation (^): Ewald's potential with conducting, or tin- 
foil, boundary conditions, Jl],|],^l]j2|] , which corresponds 
to A = in equation (pi): 
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As written this potential has a flaw: it depends on the 
arbitrary parameter a. This can be seen most easily by 
considering the mean value of over the unit cell V: 

1 Z"^*/ r^_S^ 1 / 1 ^erfc(a||x-Ln|) 
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Subtracting this mean value results in an a-independent 
potential, 
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which, furthermore, has the property that 



i J d 3 x^Jx) = Q. 



(10) 



v 



Applying the laplacian V 2 to <& Bw one verifies, after a few 
manipulations, that it is a solution of equation (0). 

Thus we can write our sought-after potential $ in the 
form 



where now 4>' satisfies the homogeneous Laplace's equa- 
tion 



In Appendix |A| the most general solution of the homoge- 
neous problem is derived. The general form found by de 
Leeuw et al., equation (Q), is recovered. In the remainder 
of the paper we choose the solution to the homogeneous 
problem equal to zero, which corresponds to setting A = 
in equation (Q). 

III. INTERACTION ENERGY OF TWO IONS. 
CONTINUUM ELECTROSTATICS APPROACH 

If two ions are placed in a periodic box filled with a 
solvent and the separation between them is large com- 
pared with typical molecular correlation lengths, the free 
energy of the system is given by the classical continuum 
electrostatic energy. The effect of the solvent is merely 
to renormalize the vacuum interaction energy (dielectric 
screening), and the complicated many body problem is 
reduced to a simpler two body problem. In this section 
we focus on this interaction energy from the point of view 
of continuum electrostatics, and in the next section we 
present a microscopic treatment of the same quantity. 
For large enough separations, of course, both treatments 
should give the same results. 

Let's therefore consider a system of two point charges, 
together with their respective neutralizing backgrounds, 
in a cubic box with periodic boundary conditions filled 
with a dielectric medium of dielectric constant e. It will 
be shown that continuum electrostatics predicts that the 
interaction energy shows a rather strong finite-size effect, 
i.e., a dependence on the linear dimension L of the box. 
In what follows we set A = so that the electric field is 
derived from a well-defined periodic potential. 

In general we can expand £ as a Fourier series 
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excluding, as discussed above, the k = term. The other 
terms can all be obtained from equation (ffl): 



ij 7T rC 



together with the Fourier space version of .E = — V0: 

27TI 

E k = -—k(f> k . (15) 

The term </>o is not fixed by equation (|7|); since it doesn't 
contribute to the physical observables (like E) this term 
is completely arbitrary. In particular, the Ewald poten- 
tial <i> Ew corresponds to <fio = 0. 

Consider now a system consisting of two charge distri- 
butions, p\ and p2, which are well separated, i.e., pj is 
zero outside of a small neighborhood of point Uj , and the 
two neighborhoods do not overlap. The energy of the sys- 
tem can be written (assuming a medium with dielectric 
constant e) 

E =^j d 3 x\\VH 2 (16) 



where <p satisfies the equation 

4-7T 

V 2 #z) = ( Pl (x) + p 2 (xj) . (17) 

£ 

For simplicity we assume that 

Qj = J d 3 x Pj ± 0, Qi + Q 2 = (18) 



otherwise a further neutralizing background must be 
added to ( |T7| ) . We want to extract from the total energy 
E the interaction piece. In an infinite system (L = oo) 
this is done by subtracting from E the terms 

E i = £ /^llv<M 2 (is) 



where 

4"7T 

VVi(aj) - -— Pj (x). (20) 

In a periodic box, as discussed previously, this equation 
must be replaced by 

V%(x) = ~^^ Pj (x)-^y (21) 

Subtracting E\ + E 2 from E we obtain the interaction 
energy, E int 

£int = ^j / '^.rV-vV.,. 122 i 



Integrating by parts and using equation ( pl| ) we can 
rewrite this as 
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Notice that 



- | d 3 2/0 2 (y) (24) 



is nothing else than the fc = term in the Fourier de- 
composition, sometimes also called the zero mode. This 
simple argument shows that to calculate the electrostatic 
interaction energy between two charge distributions the 
k = Fourier mode of the potential must be removed. 

Specializing to the case of two point charges, of charge 
+q and —q, situated at Ui and u 2 , respectively, we find 

Q 2 

E int = -— $ Ew («i -u 2 ) (25) 

£ 

where $ Ew is defined by equation (||) so that its integral 
over the unit cell vanishes identically. 



IV. INTERACTION ENERGY OF TWO IONS. 
FREE ENERGY APPROACH 

In this section we show again, but from a statistical 
mechanics perspective, that the interaction energy of two 
ions in solution coincides with the potential if and only 
if the latter satisfies that its volume average vanishes. 

Consider two ions of charges q\ and q 2 immersed in a 
molecular solvent inside a box of volume V . Under very 
general conditions it is possible to write the free energy 
of solvation of the ions as a quadratic function of the 
charges (see Appendix [§] and reference ]2q| ): 

i^solvCei) 32) = -F so iv(<7i) + F so i v (q 2 ) + qiq 2 ^(xx - x 2 ) 

+ q 1 q 2 W(x 1 -X 2 ) (26) 

where F so i v (q) is the free energy of solvation of one ion 
immersed in the solvent when the charge in the other 
is set to zero (there is still a dependence on the cavity 
produced by the other ion, but it is negligible at large 
separations of the charges). Here Xi (i = 1 and 2) are 
the positions of the ions inside the box and $ is the the 
interaction energy in the absence of solvent. The quantity 
W(x\ — x 2 ) is the solvent contribution to the effective 
interaction energy and can be written as (Appendix |b|) 

1 d 2 In Z ( 

W{x\ - x 2 ) = --- — — = -f3[ {<fr S o\v{xi)(j) sol v(x 2 )) - (^soiv(a5i))(^soiv(a52)> ) (27) 
poqidq 2 V 

where tfi so \ v (x) is the total electrostatic potential pro- 
duced by the solvent atoms at the position x. The mean 
values are computed with the ionic charges set to zero. As 
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written, W depends on the size of the ionic cavities; how- 
ever, this dependence is expected to be of short range. 
If the separation jx% — aial between the ions is larger 
than the corresponding short-range correlation length, it 
should be a good approximation to replace the mean val- 
ues in equation (|27|) by mean values computed for zero 
charge and in the absence of any cavity, which we will 
denote by the subscript 'nc' 



nc (^solv (*^1 ))nc(0solv(aJ2))»c • (28) 



By translation invariance the averages {4> S oiv(x)) nc 
are position independent and, moreover, vanish in most 
cases. We will keep them, however, for the sake of gen- 
erality If the simulation volume V were infinite we 
would expect this correlation function to decay to zero 
as ||xi — Xgfl approaches infinity. From continuum elec- 
trostatics | p4| , we expect the decay to be given by the 
asymptotic formula 

W(x 1 -x 2 )^(--l) 1 , (29) 

\e ) \x 1 - x 2 \\ 

where e is the solvent dielectric constant. Note that the 
physical decay of the correlations selects a "natural" scale 
for the interaction energy: it should be zero at infinite 
separation. In a finite, periodic box, continuum electro- 
statics leads us to expect that the Coulomb potential in 
equation (p9| ) should be replaced by its appropriate pe- 
riodic version: the Ewald potential 

W(xx-x 2 )~ Q-l^ Ew (x!-x 2 ). (30) 

In the infinite volume case we require that the potential 
go to zero at infinite separation. Since the ions in a finite 
box cannot be separated to an infinite distance, what is 
the equivalent "natural" scale for the interaction energy? 
A simple argument provides an answer, which is borne 
out by a further analysis. The interaction energy is that 
part of the energy that depends on the positions of both 
ions; therefore, any part of the potential in ( |30"| ) that does 
not depend on the positions of the ions must be removed. 
One way to do this is to set the k = Fourier compo- 
nent of the potential to zero. All the other modes depend 
in an essential way on the positions of the ions and are 
left untouched. Thus we are led to the conclusion that 
the "natural" choice for periodic boundary conditions is 
that which has no k = Fourier mode. A more rigor- 
ous way to arrive at this same conclusion is to observe 
that the integral over the simulation box vanishes (see 
Appendix [b]) : 



f 

V 



d U (0 S olv(u)0solv(O))„c - (0solv(M))nc(0solv(O))nc = 0. (31) 



This integral is in fact the k = Fourier mode of this 
two-point correlation function. 
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Since the Fourier modes depend on the simulation vol- 
ume, the requirement that the k = mode be zero in- 
troduces a subtle volume dependence. The next section 
shows with an example that this dependence can be large. 



V. INTERACTION ENERGY OF TWO IONS. 
EXPLICIT SIMULATION 

In this section we present the results of molecular dy- 
namics simulations of two ions in aqueous solution us- 
ing the Ewald summation method. After evaluating the 
charging free energy the volume dependence is discussed. 

In a large dielectric medium the contribution of the sol- 
vent to the free energy of charging, equation ([5(j|), almost 
exactly cancels the direct interaction energy (i.e., inter- 
action in the absence of solvent). In order to correctly 
model the dielectric response of the solvent in molecular 
simulations this balance needs to be accurately evaluated. 

From equation (p6|), the free energy of solvation of two 
charged particles can be related to several terms that 
are easily identifiable with particular interactions within 
the system. The first two terms represent the free en- 
ergy of solvation of each ion in solution when the other 
ion has charge zero; the second term is the direct inter- 
action contribution, and the last term is related to the 
solvent dielectric shielding. This last term is what will 
be obtained from simulations. There are at least two 
ways in which the solvent dielectric shielding contribu- 
tion can be explicitly evaluated: one is in the absence of 
the charges through a solvent correlation function (see 
equation (p27|)); the other is in presence of the charges 
by relating the dielectric shielding to the average electro- 
static potentials of the systems with two and one charged 
particles (see Appendix g). The simulations discussed in 
this section were designed to use the second approach. 

Two oppositely charged chloride-like ions were situ- 
ated on an axis perpendicular to two faces of the unit 
cell 10 A apart in a 32 A cubic box with 1103 SPC wa- 
ter molecules (the orientation is important as the Ewald 
potential is not isotropic, see Fig. |l). Three 200 picosec- 
ond molecular dynamics simulations were performed at 
room temperature for systems with the two charged par- 
ticles, and systems with only one of the particles charged. 
The average electrostatic potentials due to the solvent 
were evaluated at the cavity sites and the free energy of 
charging was calculated from the difference of the average 
potentials using the formula 

AF ^ 9i o^- 9 = ^ {(<i>soiv(x + )) q - q - (<t> so i v (x~)) q - q ) 

(32) 

where x + and x~ are, respectively, the positions of the 
positive and negative ions (see Appendix [^) . More de- 
tails about the methodology and additional studies of the 
effect of the truncation scheme on the dielectric shielding 
can be found in a related paper p9| . 
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Figure |2j shows the running average of the difference 
between the electrostatic potentials due to the solvent at 
the ionic positions. In order to obtain the free energy 
of charging due to the solvent the number obtained is 
divided by two, yielding —6.1 kcal/mol. The direct in- 
teraction energy of the ions in the absence of the solvent 
is —6.12 kcal/mol; therefore the near cancellation of the 
direct interaction by the interaction with the solvent is 
observed, as is expected for a high dielectric solvent (wa- 
ter). 

It is relevant at this point to make a comparison 
between the finite volume direct interaction energy in 
the absence of solvent, —6.12 kcal/mol, and the infi- 
nite volume direct interaction energy, —32.77 kcal/mol. 
What produces this large difference? If we compare the 
fields produced by the Coulomb and Ewald potentials 
(Fig. |3|) , a much smaller difference is found. On the other 
hand, the numbers just quoted refer to the (unmodified) 
Coulomb and the particular Ewald potential that satis- 
fies the condition of vanishing volume average (no zero 
mode). These are compared in Fig. [|. The Ewald poten- 
tial is shifted down by a considerable amount with respect 
to the Coulomb potential. Besides this shift there are 
other differences that appear because of the requirement 
of periodicity (note that the Coulomb potential does not 
meet the boundaries smoothly). These distortions are 
also of interest and will be the subject of a future publi- 
cation. 

We can estimate the importance of this shift downward 
with respect to the Coulomb potential as a function of 
the simulation volume as follows. Neglecting, to a first 
approximation, the distortions mentioned in the previ- 
ous paragraph, the amount of this shift should be given 
by the zero mode, or volume average, of the Coulomb 
potential. Actually, since the Coulomb potential is not 
periodic this computation is ill-defined; however, we can 
replace it by its simplest periodic extension, the "min- 
imum image" potential. Its volume average is easy to 
calculate and we obtain 

ill 

— / d 3 x - — r = — dx dy dz — , 

VJ ||X|| LJ J y J ^ x 2 +y 2 +z 2 

V 

~ Y x 1.19 (33) 

where L is the linear dimension of the simulation box 
(V = L 3 ). For the case considered above [L = 32 A) 
this gives about 24 kcal/mol, a substantial portion of 
the difference (32.77 — 6.12 = 26.65) between the infinite 
volume limit and the correct finite volume answer. This 
justifies, a posteriori, the neglect of the shape distortions 
as a first approximation. Thus, to a large extent, we can 
say that the Ewald potential, equation (||), is "just" the 
Coulomb potential but shifted down to guarantee that 
its mean value over the finite box is zero. This shift 
introduces a volume dependence that, as we have seen, 
is not negligible. 
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VI. DISCUSSION 



In the previous sections we have seen how a strong 
finite-size effect arises, which shifts the interaction energy 
of two point charges by a nonnegligible amount for the 
size of boxes that are customary at present. We have also 
discussed why it is mostly due to a shift downwards of 
the Ewald with respect to the Coulomb potential. 

As shown by Del Buono et al. |lj|, when this effect is 
taken into account different calculations of the dielectric 
response give results that are consistent with each other 
and with experimental observations. This is expected for 
intensive quantities, like the dielectric constant, which 
are not strongly dependent on the shape of the sample 
p6|^7t . For the calculation of extensive quantities, how- 
ever, the finite size effects described in this paper can 
produce results that are not directly comparable with ex- 
periment. A particularly interesting example is provided 
by the calculation of pK a values [ ]i1|p8| -[3l|| , or rather, the 
shift, of titratable residues in biological macromolcculcs. 
Since these shifts can be related to the interaction en- 
ergy of the set of titratable residues we expect a 
dependence on the size of the simulation box. 

The constant field in equation ([[]) introduces prob- 
lems for the consistent simulation of the system unless 
A = 0. The problems arise because of the fact that the 
system "remembers" the unit cell from which it was con- 
structed by infinite replication, and, as a consequence, 
this field does not depend continuously on the position 
of the charges: there is a finite jump whenever a charge 
crosses the boundaries of the unit cell. It is interesting 
to notice that if, instead of point charges point dipoles 
are used, these discontinuities do not arise, since the to- 
tal dipole moment is independent of the position of the 
dipoles p^l . However, in simulations using explicit atom 
models these jumps are bound to occur and it is not 
clear what the correct implementation of the dynamics 
is in that case. Until this point is clarified it seems best 



to set A = (see, however, reference 1 33 for a possible 
solution). 

The appearance of strong finite-size effects is not sur- 
prising, given the long-range nature of the electrostatic 
forces. Many earlier simulations were focused either on 
radial distribution functions, g(r), which show remark- 
able insensitivity to the detailed treatment of the elec- 
trostatic forces, or on the computation of the dielectric 
constant, which is an intrinsic measure of the response 
of the solvent and not affected by these finite-size effects. 
It is well known, however, that the Kirkwood Gx-factor, 
and its detailed relation to the dielectric constant of the 
medium, is strongly dependent on the choice of bound- 
ary conditions jj] . Nowadays many simulations, particu- 
larly those concerned with biological molecules, attempt 
to compute (differential) free energies, and it is reason- 
able to expect that these calculations will be dependent 
on the size of the simulation box through the volume de- 
pendence of the electrostatic potential. How strong this 
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dependence really is in a realistic simulation is a question 
that we will attempt to answer in the future. 

After this work was completed we received a preprint 
by Hummer, Pratt and Garcia |Q where a related finite- 
size effect was reported that affects the calculation of free 
energies of hydration of single ions. 
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APPENDIX A: 

In this appendix we derive the most general form of a 
solution to the homogeneous equations 

V E = and Vx£ = (Al) 

for the electrostatic field in a periodic box of linear di- 
mension L. 

Algebraic topology |55|,[56| can be used to prove that 
the only periodic solutions to (|Al| ) are constant fields. 
Note that if 

- V0'(x) - E , (A2) 

for a constant Eq, then <f>' must have the form 

4>'(x) = —x ■ Eo + ^-independent terms (A3) 

which is not a periodic potential. 

The explicit form of Eq given by de Leeuw et al. |2j ] 
can be recovered by the following argument. First, we 
make the natural assumption that Eq depends only on 
the distribution of charge inside the unit cell, that is, we 
assume there are no charges "at infinity." Moreover, we 
assume this dependence is linear and that the "charges 
inside the unit cell" include the homogeneous background 
(if present). Introducing the Green's function G(x, u) we 
can then write 



?o = J d 3 uG(x,u) U(w)-i J d 3 vp(v) 

V \ V J 



(A4) 



Below we will see how the form of G(x, u) is constrained 
by the requirement that Eq be constant; in the end we 
will recover the expression obtained by de Leeuw et al. 
Notice that since 
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V 



\ V 



d 6 u p{u) - — / d s v p(v) = (A5) 



for any charge distribution p, we are free to add to 
G(x, u) an arbitrary it-independent term. Furthermore, 
because of translational invariance, G can only depend 
on relative coordinates: G(x,u) = G(x — u). Taking 
derivatives with respect to on both sides of equation 



(A4), one can prove that G must be of the form 



G(x,u) = ^(x-u) + G . (A6) 

Here we have again assumed that the symmetries of the 
cubic unit cell are preserved; otherwise a slightly more 
general, albeit still linear, form for G is obtained. The 
factor of 1/V is present to make A dimensionless. The 
constant Go drops from the integral, as does the term 
proportional to x, and we are left with the final result: 

E = -^J d 3 uu L(u) - I y d 3 vp(v) J . (A7) 

V V V J 

In particular, for a collection of charges {q a } at positions 
{u a }, this expression becomes 

E Q = --^^2q a (u a -M ), (A8) 



where Uq is the center of the unit cell, given by 
1 

V 



u = i / d 3 uu. (A9) 



V 

(Notice that if Qa = Q the sum in (|A8| ) is independent 
of Uq.) Thus Eq is proportional to the dipole moment of 
the unit cell with respect to its center. The only remain- 
ing freedom lies in the choice of A which, as shown by de 
Leeuw et al., is strongly dependent on the regularization 
procedure. 

Even though Eq cannot be written as the gradient of a 
periodic function we can write Eq = — V(/> GX t, where the 
"extrinsic" potential is given by 

a 

Under the condition that ^2 a q a — 0, </> ex t can be rewrit- 
ten, up to an ^-independent term, in the form 

0«t(aO = i^^2<l a \\x-u a \ 2 . ( AU ) 



Interestingly, written in this way it takes the same form 
as the quadratic term in the expansion of the Ewald po- 
tential around the source: 



13 



2^ . „2 



<& B (x — u) = c + H la; — wll + higher-order terms. (A12) 

I a; — u\\ 3V 



This observation can be used to explain the results 
recently reported by Roberts and Schnitker j$7| and will 
be the subject of a future publication. 



It is important to realize that equation (A7) only 
makes sense after a unit cell has been chosen, which is 
then replicated ad infinitum: different choices for the unit 
cell will give different total dipole moments. Fig. |E| shows 
schematically how the choice of unit cell in a periodic 
system influences the dipole moment. The white circles 
represent positively charged atoms and the black circles 
correspond to the negative charges. As the figure shows, 
one can select different unit cells where the white circles 
are on different sides with respect to the black circles, 
therefore changing the dipole moment with respect to 
the center of the unit cell. Thus, if A ^ the system 
"remembers" the unit cell from which it was constructed 
by infinite replication. This memory of the unit cell com- 
plicates the dynamics and it is not clear to us what the 
correct treatment is. 

We should also stress that our derivation of equation 



( A7) is independent of any regularization (as long as the 
regularization respects the symmetries of the unit cell, 
see above): only the value of A can be regularization 
dependent, as shown in some specific examples by de 
Leeuw et al. [Eljl . 



APPENDIX B: 



In this Appendix we derive some identities for the 
derivatives of the free energy of a system of charges with 
respect to those charges. 

Consider a system of N fixed charges, {q a ,x a }, im- 
mersed in a solvent, constrained to a periodic box of lin- 
ear dimension L. We do not need to assume that the fixed 
charges are pointlike. Instead, the interaction with the 
solvent is taken to consist of an electrostatic piece plus 
a van der Waals-like short-range repulsion. The parti- 
tion function for the system can be decomposed into the 
product of a solvent-independent piece and the "reduced" 
partition function 



Z = Z({q a ,x a }) = Jl[du 3 exp (-pJ2 V <*,3 {x a , Uj )~f3H({u})\ (Bl) 

3 \ a,j J 



where the sum over a is over the N fixed charges and 
that over j is over the solvent atoms. The potential V a .j 
gives the interaction energy between the charge q a and 
the jth solvent atom, which we take to have the form 

V a j(x a ,Uj) = q a q.j <&{x a - Uj) + *"-(s a - uj) (B2) 

where ^> S I ' is the short-range (van der Waals) interaction 
and $ is the electrostatic potential of a unit charge. For 
the sake of notational simplicity we define 
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<j> S ow{x) = VV$(a; - uj). (B3) 



3 



The Hamiltonian contains all the energy terms that 
depend only on the solvent degrees of freedom. 

The aim of molecular mechanics simulations is to com- 
pute one or several averages with respect t o th e Boltz- 



mann probability distribution giving rise to (Bl). In this 
Appendix we will be concerned with averages that can be 
obtained as derivatives of Z with respect to the charges. 
One can also consider, of course, derivatives with respect 
to the coordinates of the fixed charges (see and ) . 

Taking a logarithmic derivative of Z with respect to 
q a we obtain 



d\nZ 

= ~P 



jY]_duj exp I -j3y^V atj (x a ,Uj) - (3H({u}) J (f> so i v (x a ) 
J l[d Uj exp —j3'Y^V a ,j(x a , Uj) — f3H{{u}) 

3 \ a,j J 



-0(<p so iy{x a )) q (B4) 



where the subscript q is meant to indicate that this mean 
value is to be computed in the presence of the N charges. 

From equation ( p34| ) and the relation In Z = —f3F be- 
tween the free energy and the partition function it follows 
that the free energy of charging one ion immersed in the 
solvent is given by 

if 

AF qi ^ qj = J dq((j) solv (x a )) q , (B5) 

with a simple extension to the case of several ions in so- 
lution. If the response to a charging process is linear 
p5| , that is, electrostriction and dielectric saturation ef- 
fects are negligible, then the mean values {4>soh/{. x a))q are 
linear in q and therefore the free energy change can be 
written simply in terms of the mean values at the end- 
points 

AF qi ^ qf = g/ - ^ (((f> S ol v {x a )) qf + (0 so lvOa)) gi ) ■ 

(B6) 

The second derivatives of the free energy have also a 
simple form: 

<pnZ = flfooiv(g/0) g = p2( {(f>sowiXa ^ solv ( x )) _ (0 solv ( a5Q )) 9 (0 solv (^)} 9 ). |B7) 
0q a dq f 3 dq a \ 1 

The higher derivatives of In Z with respect to the charges 
can be treated along the same lines but they are not 
needed in what follows. 



Note that in equation (B7) there is a hidden depen- 



dence on the positions and values of all charges, not just 
those indexed by a and /3, because of the V\j terms in 
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the Boltzmann weights. This dependence is expected to 
disappear when x a and xp are both far from all the other 
charges and from each other. In this regime the correla- 



tion function in equation (B7) should be independent of 



the value of the charges as well as the fact that they are 
not pointlike. Thus in this limit this correlation function 
can be computed by setting V\,j to zero, that is, in the 
presence of the pure solvent. Since these now are inde- 
pendent of the details of the system and only depend on 
the solvent degrees of freedom, which are integrated out, 
we can write down some general relations. 

Let's denote by (• • -) nc the mean values obtained by 
dropping these terms, i.e., 

/ Y[ duj A exp (-/3£T({u})) 

(A) nc = —pi . (B8) 

j !</".. exp (-/3ff ({«})) 



Since the exponents do not depend on the positions of 
the fixed charges, the order of operation of derivatives, 
integrals and Boltzmann averaging can be interchanged: 

— (0solv(^))nc = ( )nc (B9) 

(JJL j (J JL j 



d 3 X (0 S olvO))nc = ( / d 3 x0 solv (x)) nc . (BIO) 
In the same way we obtain the relation 

d 3 X ^(0solv(O)(/) S ol v (a;))nc - (0solv(O))nc(0solv(a;))i 

(^soiv(O) J d 3 x(t) soW (x)) DC - (0 so iv(O)) nc ( / d a x^ l]v (x)) lh . (Bll ) 

Since 



J d 3 x 4> so i v (x) =y^ q j j d 3 x $(x - Uj) 



J2lA / d 3 x<P(x) (B12) 



3 



V 



is independent of the position of the solvent atoms it 
follows that 



WO) / d 6 x<p solv (x)) ac = (B13) 

(0)) nc ( / d 3 X 0solv(*))nc- 



V 

From this it follows that 



d x (0 so iv(O)0 so iv(a;)) nc (^solv (O))„ c (0soiv(a;))„c =0. (B14) 

The mean values 



16 



(</> so 1v(:ei)<?Wv(:e2) ■ ■ ■ 0soiv(a5;v))nc (B15) 

are independent of any details of the fixed charges, since 
they are computed in the presence of the pure solvent. 
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FIG. 1. Contour plot of the Ewald potential. The labels correspond to the levels in kcal/mol-au, for a simulation box of 
32 A. The contours lie all in the plane z = 0. 

FIG. 2. Running average of the difference in the electrostatic potential produced by the solvent in the charged and uncharged 
states. 

FIG. 3. Comparison between Ewald (full curve) and Coulomb (dashed curve) fields for a unit charge at the origin in a box 
of side length 32 A. 



FIG. 4. Comparison between Ewald (full curve) and Coulomb (dashed curve) potentials for a unit charge at the origin in a 
box of side length 32 A. 



FIG. 5. Two choices of unit cell in a periodic system. 
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